canvod.grids API Reference¶

class BaseGridBuilder(ABC):
    """Abstract base for hemispherical grid builders.

    Parameters
    ----------
    angular_resolution : float
        Angular resolution in degrees
    cutoff_theta : float
        Maximum polar angle cutoff in degrees
    phi_rotation : float
        Rotation angle in degrees (applied to all phi values)

    """

    def __init__(
        self,
        angular_resolution: float = 2,
        cutoff_theta: float = 0,
        phi_rotation: float = 0,
    ) -> None:
        """Initialize the grid builder.

        Parameters
        ----------
        angular_resolution : float, default 2
            Angular resolution in degrees.
        cutoff_theta : float, default 0
            Maximum polar angle cutoff in degrees.
        phi_rotation : float, default 0
            Rotation angle in degrees.

        """
        self.angular_resolution = angular_resolution
        self.angular_resolution_rad = np.deg2rad(angular_resolution)
        self.cutoff_theta = cutoff_theta
        self.cutoff_theta_rad = np.deg2rad(cutoff_theta)
        self.phi_rotation = phi_rotation
        self.phi_rotation_rad = np.deg2rad(phi_rotation)
        self._logger = _get_logger()

    @abstractmethod
    def _build_grid(
        self,
    ) -> (
        tuple[pl.DataFrame, np.ndarray, list[np.ndarray], list[np.ndarray]]
        | tuple[
            pl.DataFrame, np.ndarray, list[np.ndarray], list[np.ndarray], dict[str, Any]
        ]
    ):
        """Build grid.

        Returns
        -------
        grid : pl.DataFrame
            Grid cells
        theta_lims : np.ndarray
            Theta band limits
        phi_lims : list[np.ndarray]
            Phi limits per band
        cell_ids : list[np.ndarray]
            Cell IDs per band
        extra_kwargs : dict, optional
            Additional metadata

        """

    @abstractmethod
    def get_grid_type(self) -> str:
        """Get grid type identifier."""

    def build(self) -> GridData:
        """Build hemisphere grid.

        Returns
        -------
        GridData
            Complete grid data structure

        """
        self._logger.info(
            "grid_build_started",
            grid_type=self.get_grid_type(),
            angular_resolution=self.angular_resolution,
        )

        result = self._build_grid()

        if len(result) == 4:
            grid, theta_lims, phi_lims, cell_ids = result
            extra_kwargs = {}
        elif len(result) == 5:
            grid, theta_lims, phi_lims, cell_ids, extra_kwargs = result
        else:
            raise ValueError(f"Invalid grid builder result: {len(result)} elements")

        # Apply phi rotation if specified (vectorized operations)
        if self.phi_rotation_rad != 0:
            grid = grid.with_columns(
                [(pl.col("phi") + self.phi_rotation_rad) % (2 * np.pi)]
            )

            if "phi_min" in grid.columns:
                grid = grid.with_columns(
                    [
                        (
                            (pl.col("phi_min") + self.phi_rotation_rad) % (2 * np.pi)
                        ).alias("phi_min"),
                        (
                            (pl.col("phi_max") + self.phi_rotation_rad) % (2 * np.pi)
                        ).alias("phi_max"),
                    ]
                )

        self._logger.info("grid_build_complete", ncells=len(grid))

        return GridData(
            grid=grid,
            theta_lims=theta_lims,
            phi_lims=phi_lims,
            cell_ids=cell_ids,
            grid_type=self.get_grid_type(),
            **extra_kwargs,
        )

@dataclass(frozen=True)
class GridData:
    """Immutable container for hemispherical grid structure.

    Parameters
    ----------
    grid : pl.DataFrame
        Grid cells with phi, theta, and bounds
    theta_lims : np.ndarray
        Theta band limits
    phi_lims : list[np.ndarray]
        Phi limits per theta band
    cell_ids : list[np.ndarray]
        Cell IDs per theta band
    grid_type : str
        Grid type identifier
    solid_angles : np.ndarray, optional
        Solid angles per cell [steradians]
    metadata : dict, optional
        Additional grid metadata
    voronoi : Any, optional
        Voronoi tessellation object (for Fibonacci grids)
    vertices : np.ndarray, optional
        3D vertices (for triangular grids)
    points_xyz : np.ndarray, optional
        3D point cloud (for Fibonacci grids)
    vertex_phi : np.ndarray, optional
        Vertex phi coordinates
    vertex_theta : np.ndarray, optional
        Vertex theta coordinates

    """

    grid: pl.DataFrame
    theta_lims: np.ndarray
    phi_lims: list[np.ndarray]
    cell_ids: list[np.ndarray]
    grid_type: str
    solid_angles: np.ndarray | None = None
    metadata: dict | None = None
    voronoi: Any | None = None
    vertices: np.ndarray | None = None
    points_xyz: np.ndarray | None = None
    vertex_phi: np.ndarray | None = None
    vertex_theta: np.ndarray | None = None

    @property
    def coords(self) -> pl.DataFrame:
        """Get cell coordinates."""
        return self.grid.select(["phi", "theta"])

    @property
    def ncells(self) -> int:
        """Number of cells in grid."""
        return len(self.grid)

    def get_patches(self) -> pl.Series:
        """Create matplotlib patches for polar visualization."""
        patches = [
            Rectangle(
                (row["phi_min"], row["theta_min"]),
                row["phi_max"] - row["phi_min"],
                row["theta_max"] - row["theta_min"],
                fill=True,
            )
            for row in self.grid.iter_rows(named=True)
        ]
        return pl.Series("Patches", patches)

    def get_solid_angles(self) -> np.ndarray:
        """Calculate solid angle for each cell [steradians]."""
        if self.solid_angles is not None:
            return self.solid_angles

        # HEALPix
        if self.grid_type == "healpix" and "healpix_nside" in self.grid.columns:
            try:
                import healpy as hp

                nside = int(self.grid["healpix_nside"][0])
                return np.full(
                    len(self.grid), hp.nside2pixarea(nside), dtype=np.float64
                )
            except ImportError:
                pass

        # Geodesic
        if self.grid_type == "geodesic" and "geodesic_vertices" in self.grid.columns:
            return self._compute_geodesic_solid_angles()

        # HTM
        if self.grid_type == "htm" and "htm_vertex_0" in self.grid.columns:
            return self._compute_htm_solid_angles()

        # Fibonacci
        if self.grid_type == "fibonacci" and "voronoi_region" in self.grid.columns:
            return self._compute_voronoi_solid_angles()

        # Default
        return self._geometric_solid_angles()

    def _compute_htm_solid_angles(self) -> np.ndarray:
        """Compute solid angles for HTM triangular cells."""
        solid_angles = []

        for row in self.grid.iter_rows(named=True):
            v0 = np.array(row["htm_vertex_0"])
            v1 = np.array(row["htm_vertex_1"])
            v2 = np.array(row["htm_vertex_2"])

            # Spherical excess formula
            a = np.arccos(np.clip(np.dot(v1, v2), -1, 1))
            b = np.arccos(np.clip(np.dot(v0, v2), -1, 1))
            c = np.arccos(np.clip(np.dot(v0, v1), -1, 1))

            s = (a + b + c) / 2
            tan_E_4 = np.sqrt(
                np.tan(s / 2)
                * np.tan((s - a) / 2)
                * np.tan((s - b) / 2)
                * np.tan((s - c) / 2)
            )
            E = 4 * np.arctan(tan_E_4)

            solid_angles.append(E)

        return np.array(solid_angles)

    def _compute_geodesic_solid_angles(self) -> np.ndarray:
        """Compute solid angles for geodesic triangular cells."""
        vertices = self.vertices
        if vertices is None:
            return self._geometric_solid_angles()

        solid_angles = []
        for row in self.grid.iter_rows(named=True):
            v_indices = row["geodesic_vertices"]
            v0, v1, v2 = vertices[v_indices]

            a = np.arccos(np.clip(np.dot(v1, v2), -1, 1))
            b = np.arccos(np.clip(np.dot(v0, v2), -1, 1))
            c = np.arccos(np.clip(np.dot(v0, v1), -1, 1))

            s = (a + b + c) / 2
            tan_E_4 = np.sqrt(
                np.tan(s / 2)
                * np.tan((s - a) / 2)
                * np.tan((s - b) / 2)
                * np.tan((s - c) / 2)
            )
            E = 4 * np.arctan(tan_E_4)

            solid_angles.append(E)

        return np.array(solid_angles)

    def _compute_voronoi_solid_angles(self) -> np.ndarray:
        """Compute solid angles for Voronoi cells."""
        if self.voronoi is None:
            return self._geometric_solid_angles()

        sv = self.voronoi
        solid_angles = []
        for row in self.grid.iter_rows(named=True):
            region = row["voronoi_region"]
            if len(region) < 3:
                solid_angles.append(np.nan)
                continue

            vertices = sv.vertices[region]
            center = np.array(
                [
                    np.sin(row["theta"]) * np.cos(row["phi"]),
                    np.sin(row["theta"]) * np.sin(row["phi"]),
                    np.cos(row["theta"]),
                ]
            )

            total_angle = 0
            n = len(vertices)
            for i in range(n):
                v1 = vertices[i]
                v2 = vertices[(i + 1) % n]
                a = np.arccos(np.clip(np.dot(center, v1), -1, 1))
                b = np.arccos(np.clip(np.dot(center, v2), -1, 1))
                c = np.arccos(np.clip(np.dot(v1, v2), -1, 1))
                s = (a + b + c) / 2
                tan_E_4 = np.sqrt(
                    np.tan(s / 2)
                    * np.tan((s - a) / 2)
                    * np.tan((s - b) / 2)
                    * np.tan((s - c) / 2)
                )
                E = 4 * np.arctan(tan_E_4)
                total_angle += E
            solid_angles.append(total_angle)

        return np.array(solid_angles)

    def _geometric_solid_angles(self) -> np.ndarray:
        """Fallback geometric calculation."""
        solid_angles = []
        for row in self.grid.iter_rows(named=True):
            delta_phi = row["phi_max"] - row["phi_min"]
            cos_diff = np.cos(row["theta_min"]) - np.cos(row["theta_max"])
            omega = delta_phi * cos_diff
            solid_angles.append(omega)
        return np.array(solid_angles)

    def get_grid_stats(self) -> dict:
        """Get grid statistics including solid angle uniformity."""
        solid_angles = self.get_solid_angles()

        stats = {
            "total_cells": self.ncells,
            "grid_type": self.grid_type,
            "theta_bands": len(self.theta_lims),
            "cells_per_band": [len(ids) for ids in self.cell_ids],
            "solid_angle_mean_sr": float(np.mean(solid_angles)),
            "solid_angle_std_sr": float(np.std(solid_angles)),
            "solid_angle_cv_percent": float(
                np.std(solid_angles) / np.mean(solid_angles) * 100
            ),
            "total_solid_angle_sr": float(np.sum(solid_angles)),
            "hemisphere_solid_angle_sr": 2 * np.pi,
        }

        # Add HEALPix-specific info
        if self.grid_type == "healpix" and "healpix_nside" in self.grid.columns:
            try:
                import healpy as hp

                nside = int(self.grid["healpix_nside"][0])
                stats["healpix_nside"] = nside
                stats["healpix_npix_total"] = hp.nside2npix(nside)
                stats["healpix_pixel_area_sr"] = hp.nside2pixarea(nside)
                stats["healpix_resolution_arcmin"] = hp.nside2resol(nside, arcmin=True)
            except ImportError:
                pass

        return stats

class GridType(Enum):
    """Available grid projection types for hemispherical tessellation."""

    EQUAL_AREA = "equal_area"  # Equal solid angle (ring-based)
    EQUAL_ANGLE = "equal_angle"  # Equal angular spacing
    EQUIRECTANGULAR = "equirectangular"  # Simple rectangular
    HEALPIX = "healpix"  # Hierarchical equal area
    GEODESIC = "geodesic"  # Icosahedral triangular
    FIBONACCI = "fibonacci"  # Golden spiral + Voronoi
    HTM = "htm"  # Hierarchical Triangular Mesh

class EqualAreaBuilder(BaseGridBuilder):
    """Equal solid angle tessellation using concentric theta bands.

    The hemisphere is divided into annular bands of constant width in theta.
    Within each band the number of azimuthal (phi) sectors is chosen so that
    every cell subtends approximately the same solid angle.  This is the only
    grid type that has been validated for scientific use in this codebase.

    Coordinate convention (physics / GNSS)
    ---------------------------------------
    * phi  ∈ [0, 2π)  – azimuthal angle from North, clockwise (navigation convention)
    * theta ∈ [0, π/2] – polar angle measured from zenith (0 = straight up,
      π/2 = horizon)

    What ``angular_resolution`` means
    ----------------------------------
    ``angular_resolution`` (degrees) sets the **width of each theta band**.
    All bands have this same width Δθ.  The *azimuthal* width of cells varies
    by band: near the zenith cells are wide in phi; near the horizon they are
    narrow, so that the solid angle stays constant.

    Mathematical construction
    -------------------------
    1. **Target solid angle** per cell is chosen equal to the solid angle of a
       cap of half-angle Δθ/2::

           Ω_target = 2π (1 − cos(Δθ/2))

    2. **Zenith cap** – a single cell covers [0, Δθ/2] in theta and the full
       azimuth [0, 2π).

    3. **Theta bands** – edges are placed at Δθ/2, 3Δθ/2, 5Δθ/2, … up to
       π/2 − cutoff_theta.  For each band [θ_inner, θ_outer] the band's
       total solid angle is::

           Ω_band = 2π (cos θ_inner − cos θ_outer)

    4. **Phi divisions** – the number of sectors in the band is::

           n_phi = round(Ω_band / Ω_target)

       Each sector spans Δφ = 2π / n_phi.  The cell centre is placed at the
       geometric midpoint of its (phi, theta) rectangle.

    Parameters
    ----------
    angular_resolution : float
        Theta-band width in degrees.  Controls both the radial resolution and
        (indirectly, via the equal-area constraint) the azimuthal resolution.
    cutoff_theta : float
        Minimum elevation above the horizon in degrees.  Bands whose outer
        edge is at or below this cutoff are omitted.  In GNSS terms this is
        the satellite elevation mask angle.
    phi_rotation : float
        Rigid rotation applied to all phi coordinates after grid construction,
        in degrees.

    """

    def get_grid_type(self) -> str:
        """Return the grid-type identifier string.

        Returns
        -------
        str
            ``"equal_area"``

        """
        return GridType.EQUAL_AREA.value

    def _build_grid(
        self,
    ) -> tuple[pl.DataFrame, np.ndarray, list[np.ndarray], list[np.ndarray]]:
        """Construct the equal-area hemisphere grid.

        Returns
        -------
        grid : pl.DataFrame
            One row per cell with columns: phi, theta, phi_min, phi_max,
            theta_min, theta_max, cell_id.
        theta_lims : np.ndarray
            Outer theta edge of each band (radians).
        phi_lims : list[np.ndarray]
            Array of phi_min values for each band.
        cell_ids : list[np.ndarray]
            Cell-id arrays, one per band.

        """
        # Theta band edges (from zenith to horizon)
        max_theta = np.pi / 2  # horizon
        theta_edges = np.arange(
            self.angular_resolution_rad / 2,
            max_theta - self.cutoff_theta_rad,
            self.angular_resolution_rad,
        )

        # Target solid angle per cell
        target_omega = 2 * np.pi * (1 - np.cos(self.angular_resolution_rad / 2))

        cells = []
        theta_lims = []
        phi_lims = []
        cell_ids = []

        # Zenith cell (special case) - only if cutoff allows
        next_cell_id = 0
        zenith_theta_max = self.angular_resolution_rad / 2

        if self.cutoff_theta_rad < zenith_theta_max:
            cells.append(
                pl.DataFrame(
                    {
                        "phi": [0.0],
                        "theta": [0.0],
                        "phi_min": [0.0],
                        "phi_max": [2 * np.pi],
                        "theta_min": [max(0.0, self.cutoff_theta_rad)],
                        "theta_max": [zenith_theta_max],
                    }
                )
            )
            theta_lims.append(zenith_theta_max)
            phi_lims.append(np.array([0.0]))
            cell_ids.append(np.array([0]))
            next_cell_id = 1

        # Build theta bands
        for iband, theta_outer in enumerate(theta_edges[1:]):
            theta_inner = theta_edges[iband]

            # Skip bands below cutoff
            if theta_outer <= self.cutoff_theta_rad:
                continue

            # Solid angle of this band
            band_omega = 2 * np.pi * (np.cos(theta_inner) - np.cos(theta_outer))

            # Number of phi divisions
            n_phi = max(1, round(band_omega / target_omega))
            phi_span = 2 * np.pi / n_phi

            cell_id_list = list(range(next_cell_id, next_cell_id + n_phi))
            next_cell_id = cell_id_list[-1] + 1

            # Use arange for better precision than linspace
            phi_min_arr = np.arange(n_phi) * phi_span
            phi_max_arr = (np.arange(n_phi) + 1) * phi_span
            phi_max_arr[-1] = 2 * np.pi  # Force exact closure

            cells.append(
                pl.DataFrame(
                    {
                        "phi": (phi_min_arr + phi_max_arr) / 2,
                        "theta": np.full(n_phi, (theta_inner + theta_outer) / 2),
                        "phi_min": phi_min_arr,
                        "phi_max": phi_max_arr,
                        "theta_min": np.full(n_phi, theta_inner),
                        "theta_max": np.full(n_phi, theta_outer),
                    }
                )
            )

            theta_lims.append(theta_outer)
            phi_lims.append(phi_min_arr)
            cell_ids.append(np.array(cell_id_list))

        if len(cells) == 0:
            raise ValueError(
                "No cells generated - check cutoff_theta and angular_resolution"
            )

        grid = pl.concat(cells).with_columns(pl.int_range(0, pl.len()).alias("cell_id"))

        return grid, np.array(theta_lims), phi_lims, cell_ids

class EqualAngleBuilder(BaseGridBuilder):
    """Equal angular spacing in both theta and phi (NOT equal area).

    Every cell is a rectangle of the same angular size Δθ × Δφ in the
    (theta, phi) parameter space.  Because solid angle depends on cos(theta),
    cells near the zenith subtend *more* solid angle than cells near the
    horizon.  This makes the grid biased toward the zenith for any
    solid-angle-weighted statistic.  **Not recommended for scientific
    analysis** – use ``EqualAreaBuilder`` instead.

    Coordinate convention (physics / GNSS)
    ---------------------------------------
    * phi  ∈ [0, 2π)  – azimuthal angle from North, clockwise (navigation convention)
    * theta ∈ [0, π/2] – polar angle from zenith

    What ``angular_resolution`` means
    ----------------------------------
    ``angular_resolution`` (degrees) is used as **both** the theta-band width
    and the phi-sector width.  The number of phi divisions is constant across
    all bands::

        n_phi = round(2π / Δθ)

    and does not change with latitude.

    Mathematical construction
    -------------------------
    1. A zenith cap cell covers [0, Δθ/2] × [0, 2π).
    2. Theta band edges are placed at Δθ/2, 3Δθ/2, … up to π/2.
    3. Within every band, the full azimuth is split into ``n_phi`` sectors of
       equal width Δφ = 2π / n_phi.
    4. Cell centres are at the midpoint of each (phi, theta) rectangle.

    Parameters
    ----------
    angular_resolution : float
        Angular spacing in degrees, applied identically in both theta and phi.
    cutoff_theta : float
        Elevation mask angle in degrees (bands below this are omitted).
    phi_rotation : float
        Rigid azimuthal rotation applied after construction, in degrees.

    """

    def get_grid_type(self) -> str:
        """Return the grid-type identifier string.

        Returns
        -------
        str
            ``"equal_angle"``

        """
        return GridType.EQUAL_ANGLE.value

    def _build_grid(
        self,
    ) -> tuple[pl.DataFrame, np.ndarray, list[np.ndarray], list[np.ndarray]]:
        """Construct the equal-angle hemisphere grid.

        Returns
        -------
        grid : pl.DataFrame
            One row per cell.
        theta_lims : np.ndarray
            Outer theta edge of each band (radians).
        phi_lims : list[np.ndarray]
            Array of phi_min values for each band (identical across bands).
        cell_ids : list[np.ndarray]
            Cell-id arrays, one per band.

        """
        max_theta = np.pi / 2
        theta_edges = np.arange(
            self.angular_resolution_rad / 2,
            max_theta - self.cutoff_theta_rad,
            self.angular_resolution_rad,
        )

        n_phi_divisions = int(2 * np.pi / self.angular_resolution_rad)

        cells = []
        theta_lims = []
        phi_lims = []
        cell_ids = []

        # Zenith
        cells.append(
            pl.DataFrame(
                {
                    "phi": [0.0],
                    "theta": [0.0],
                    "phi_min": [0.0],
                    "phi_max": [2 * np.pi],
                    "theta_min": [0.0],
                    "theta_max": [self.angular_resolution_rad / 2],
                }
            )
        )
        theta_lims.append(self.angular_resolution_rad / 2)
        phi_lims.append(np.array([0.0]))
        cell_ids.append(np.array([0]))
        next_cell_id = 1

        for iband, theta_outer in enumerate(theta_edges[1:]):
            theta_inner = theta_edges[iband]
            phi_span = 2 * np.pi / n_phi_divisions

            cell_id_list = list(range(next_cell_id, next_cell_id + n_phi_divisions))
            next_cell_id = cell_id_list[-1] + 1

            phi_min_arr = np.linspace(0, 2 * np.pi - phi_span, n_phi_divisions)
            phi_max_arr = np.concatenate((phi_min_arr[1:], [2 * np.pi]))

            cells.append(
                pl.DataFrame(
                    {
                        "phi": (phi_min_arr + phi_max_arr) / 2,
                        "theta": np.full(
                            n_phi_divisions,
                            (theta_inner + theta_outer) / 2,
                        ),
                        "phi_min": phi_min_arr,
                        "phi_max": phi_max_arr,
                        "theta_min": np.full(n_phi_divisions, theta_inner),
                        "theta_max": np.full(n_phi_divisions, theta_outer),
                    }
                )
            )

            theta_lims.append(theta_outer)
            phi_lims.append(phi_min_arr)
            cell_ids.append(np.array(cell_id_list))

        grid = pl.concat(cells).with_row_index("cell_id")
        return grid, np.array(theta_lims), phi_lims, cell_ids

class EquirectangularBuilder(BaseGridBuilder):
    """Simple rectangular grid in (theta, phi) space.

    The hemisphere is divided into a regular rectangular array: a constant
    number of theta bands, each containing the same constant number of phi
    sectors.  Every cell is an identical rectangle in angular coordinates.
    This is *structurally* identical to ``EqualAngleBuilder`` except for one
    difference in the zenith treatment: ``EqualAngleBuilder`` collapses the
    first band into a single zenith cap, while this builder does not — every
    band has the same number of sectors.

    Because solid angle depends on cos(theta), cells near the zenith subtend
    *more* solid angle than cells near the horizon.  This makes the grid
    biased toward the zenith for any solid-angle-weighted statistic.
    **Not recommended for scientific analysis** – use ``EqualAreaBuilder``
    instead.

    Coordinate convention (physics / GNSS)
    ---------------------------------------
    * phi  ∈ [0, 2π)  – azimuthal angle from North, clockwise (navigation convention)
    * theta ∈ [0, π/2] – polar angle from zenith (0 = straight up,
      π/2 = horizon)

    What ``angular_resolution`` means
    ----------------------------------
    ``angular_resolution`` (degrees) is used as **both** the theta-band width
    *and* the phi-sector width.  The grid is therefore square in angular
    coordinates::

        n_theta = round((π/2 − cutoff) / Δθ)
        n_phi   = round(2π / Δθ)
        total cells = n_theta × n_phi

    Mathematical construction
    -------------------------
    1. Theta edges are placed at ``cutoff_theta``, ``cutoff_theta + Δθ``,
       ``cutoff_theta + 2Δθ``, … up to π/2.
    2. Phi edges are placed at 0, Δθ, 2Δθ, … up to 2π.
    3. Every (theta_band, phi_sector) combination produces one cell.  The
       cell centre is the midpoint of the rectangle.
    4. No special zenith cap is created; the band nearest the zenith has
       the same number of phi sectors as all other bands.

    Parameters
    ----------
    angular_resolution : float
        Angular spacing in degrees, applied identically in both theta and phi.
    cutoff_theta : float
        Elevation mask angle in degrees.  Bands whose *inner* edge is at or
        below ``π/2 − cutoff_theta`` are omitted.
    phi_rotation : float
        Rigid azimuthal rotation applied after construction, in degrees.

    """

    def get_grid_type(self) -> str:
        """Return the grid-type identifier string.

        Returns
        -------
        str
            ``"equirectangular"``

        """
        return GridType.EQUIRECTANGULAR.value

    def _build_grid(
        self,
    ) -> tuple[pl.DataFrame, np.ndarray, list[np.ndarray], list[np.ndarray]]:
        """Construct the equirectangular hemisphere grid.

        Returns
        -------
        grid : pl.DataFrame
            One row per cell with columns: phi, theta, phi_min, phi_max,
            theta_min, theta_max, cell_id.
        theta_lims : np.ndarray
            Inner theta edge of each band (radians).
        phi_lims : list[np.ndarray]
            Array of phi_min values for each band (identical across bands).
        cell_ids : list[np.ndarray]
            Cell-id arrays, one per band.

        """
        max_theta = np.pi / 2

        theta_edges = np.arange(
            self.cutoff_theta_rad,
            max_theta + self.angular_resolution_rad,
            self.angular_resolution_rad,
        )
        phi_edges = np.arange(
            0, 2 * np.pi + self.angular_resolution_rad, self.angular_resolution_rad
        )

        cells = []

        for i in range(len(theta_edges) - 1):
            theta_min, theta_max = theta_edges[i], theta_edges[i + 1]

            for j in range(len(phi_edges) - 1):
                phi_min, phi_max = phi_edges[j], phi_edges[j + 1]

                cells.append(
                    {
                        "phi": (phi_min + phi_max) / 2,
                        "theta": (theta_min + theta_max) / 2,
                        "phi_min": phi_min,
                        "phi_max": min(2 * np.pi, phi_max),
                        "theta_min": theta_min,
                        "theta_max": theta_max,
                    }
                )

        grid = pl.DataFrame(cells).with_columns(
            pl.int_range(0, pl.len()).alias("cell_id")
        )

        theta_lims = theta_edges[:-1]
        phi_lims = [phi_edges[:-1] for _ in range(len(theta_edges) - 1)]
        cell_ids_list = [
            np.arange(i * (len(phi_edges) - 1), (i + 1) * (len(phi_edges) - 1))
            for i in range(len(theta_edges) - 1)
        ]

        return grid, theta_lims, phi_lims, cell_ids_list

class HEALPixBuilder(BaseGridBuilder):
    """HEALPix tessellation (Hierarchical Equal Area isoLatitude Pixelization).

    HEALPix partitions the sphere into 12 base pixels arranged at equal
    latitudes.  Each base pixel is recursively subdivided into 4 children,
    producing ``12 × nside²`` pixels on the full sphere, all with *exactly*
    the same solid angle.  This strict equal-area property makes HEALPix
    the gold standard for pixelisations that must be unbiased under
    solid-angle weighting.

    This builder delegates the pixel geometry entirely to the ``healpy``
    library.  It filters the full-sphere pixelisation down to the northern
    hemisphere and stores approximate bounding boxes (``phi_min/max``,
    ``theta_min/max``) derived from the pixel resolution.  The bounding
    boxes are **not** the true pixel boundaries (which are curvilinear);
    they are only approximations suitable for quick spatial queries.  For
    exact pixel membership use ``healpy.ang2pix`` directly.

    Coordinate convention
    ---------------------
    HEALPix natively uses colatitude ``theta ∈ [0, π]`` (0 = North Pole)
    and longitude ``phi ∈ [0, 2π)``.  This matches the GNSS convention used
    elsewhere in canvodpy: theta = 0 is the zenith, theta = π/2 is the
    horizon.  **No coordinate transform is applied.**

    What ``nside`` (resolution) means
    ----------------------------------
    ``nside`` is the single resolution parameter of HEALPix.  It must be a
    power of 2.  The key derived quantities are::

        n_pixels   = 12 × nside²           (full sphere)
        pixel_area = 4π / n_pixels          (steradians, exact)
        resolution ≈ √(pixel_area)          (approximate angular diameter)
                   ≈ 58.6° / nside         (degrees)

    | nside | Pixels (full) | Approx resolution | Pixel area (sr) |
    |-------|---------------|-------------------|-----------------|
    | 1     | 12            | 58.6°             | 1.049           |
    | 2     | 48            | 29.3°             | 0.262           |
    | 4     | 192           | 14.7°             | 0.065           |
    | 8     | 768           | 7.3°              | 0.016           |
    | 16    | 3 072         | 3.7°              | 0.004           |
    | 32    | 12 288        | 1.8°              | 0.001           |

    When ``nside`` is not provided, it is estimated from ``angular_resolution``
    and rounded to the nearest power of 2::

        nside_estimate = round_to_pow2( √(3/π) × 60 / angular_resolution )

    Mathematical construction
    -------------------------
    HEALPix construction is performed entirely by ``healpy``.  At a high
    level:

    1. The sphere is divided into 12 congruent base pixels (a curvilinear
       quadrilateral arrangement at three latitude zones: polar caps and
       equatorial belt).
    2. Each base pixel is subdivided into ``nside²`` equal-area children
       using a hierarchical quadtree.
    3. Pixel centres are returned by ``healpy.pix2ang(nside, ipix)`` in
       RING ordering (pixels ordered by increasing colatitude).
    4. This builder keeps only pixels with ``theta ≤ π/2 − cutoff_theta``
       (northern hemisphere above the elevation mask).

    Parameters
    ----------
    angular_resolution : float
        Approximate angular resolution in degrees.  Used only to derive
        ``nside`` when that parameter is not given explicitly.
    cutoff_theta : float
        Elevation mask angle in degrees.  Pixels with colatitude
        ``theta > π/2 − cutoff_theta`` (i.e. below the mask) are excluded.
    nside : int or None
        HEALPix resolution parameter.  Must be a power of 2.  If ``None``,
        estimated from ``angular_resolution``.
    phi_rotation : float
        Rigid azimuthal rotation applied after construction, in degrees.

    Raises
    ------
    ImportError
        If ``healpy`` is not installed.
    ValueError
        If ``nside`` is not a power of 2.

    """

    def __init__(
        self,
        angular_resolution: float = 2,
        cutoff_theta: float = 0,
        nside: int | None = None,
        phi_rotation: float = 0,
    ) -> None:
        """Initialize the HEALPix grid builder.

        Parameters
        ----------
        angular_resolution : float, default 2
            Angular resolution in degrees.
        cutoff_theta : float, default 0
            Maximum polar angle cutoff in degrees.
        nside : int | None, optional
            HEALPix nside parameter.
        phi_rotation : float, default 0
            Rotation angle in degrees.

        """
        super().__init__(angular_resolution, cutoff_theta, phi_rotation)

        # Determine nside
        if nside is None:
            nside_estimate = int(np.sqrt(3 / np.pi) * 60 / angular_resolution)
            self.nside = 2 ** max(0, int(np.round(np.log2(nside_estimate))))
        else:
            if nside < 1 or (nside & (nside - 1)) != 0:
                raise ValueError(f"nside must be a power of 2, got {nside}")
            self.nside = nside

        # Import healpy
        try:
            import healpy as hp

            self.hp = hp
        except ImportError:
            raise ImportError(
                "healpy is required for HEALPix grid. Install with: pip install healpy"
            )

        pixel_size_arcmin = self.hp.nside2resol(self.nside, arcmin=True)
        self.actual_angular_resolution = pixel_size_arcmin / 60.0

        self._logger.info(
            f"HEALPix: nside={self.nside}, "
            f"requested_res={angular_resolution:.2f}°, "
            f"actual_res={self.actual_angular_resolution:.2f}°"
        )

    def get_grid_type(self) -> str:
        """Return the grid-type identifier string.

        Returns
        -------
        str
            ``"healpix"``

        """
        return GridType.HEALPIX.value

    def _build_grid(
        self,
    ) -> tuple[pl.DataFrame, np.ndarray, list[np.ndarray], list[np.ndarray]]:
        """Build HEALPix grid for the northern hemisphere.

        Iterates over all ``12 × nside²`` pixels, retains those with
        ``theta ≤ π/2 − cutoff_theta``, and constructs approximate
        bounding boxes from the pixel resolution.

        Returns
        -------
        grid : pl.DataFrame
            One row per pixel.  Contains phi, theta (centre), approximate
            bounding-box limits, ``healpix_ipix`` (RING-ordered pixel index),
            and ``healpix_nside``.
        theta_lims : np.ndarray
            Synthetic evenly-spaced theta limits (interface compatibility only).
        phi_lims : list[np.ndarray]
            Synthetic evenly-spaced phi limits (interface compatibility only).
        cell_ids : list[np.ndarray]
            Single-element list containing the valid pixel indices.

        """
        npix = self.hp.nside2npix(self.nside)

        cells = []
        valid_pixels = []

        for ipix in range(npix):
            theta, phi = self.hp.pix2ang(self.nside, ipix)

            # Keep only northern hemisphere above the elevation mask
            if theta > (np.pi / 2 - self.cutoff_theta_rad):
                continue

            pixel_radius = self.hp.nside2resol(self.nside)

            cells.append(
                {
                    "phi": float(phi),
                    "theta": float(theta),
                    "phi_min": float(max(0, phi - pixel_radius / 2)),
                    "phi_max": float(min(2 * np.pi, phi + pixel_radius / 2)),
                    "theta_min": float(max(0, theta - pixel_radius / 2)),
                    "theta_max": float(min(np.pi / 2, theta + pixel_radius / 2)),
                    "healpix_ipix": int(ipix),
                    "healpix_nside": int(self.nside),
                }
            )
            valid_pixels.append(int(ipix))

        if len(cells) == 0:
            raise ValueError("No valid HEALPix pixels found in hemisphere")

        grid = pl.DataFrame(cells)

        grid = grid.with_columns(
            [
                pl.col("healpix_ipix").cast(pl.Int64),
                pl.col("healpix_nside").cast(pl.Int64),
            ]
        )

        theta_unique = sorted(grid["theta"].unique())
        n_theta_bands = len(theta_unique)

        # NOTE: These limits are SYNTHETIC and do NOT correspond to actual
        # HEALPix pixel boundaries. They exist only for interface
        # compatibility with ring-based grids. For spatial queries, use the
        # per-pixel theta_min/max and phi_min/max columns instead.
        theta_lims = np.linspace(0, np.pi / 2, min(n_theta_bands, 20))
        phi_lims = [np.linspace(0, 2 * np.pi, 20) for _ in range(len(theta_lims))]

        cell_ids_list = [np.array(valid_pixels, dtype=np.int64)]

        return grid, theta_lims, phi_lims, cell_ids_list

    def get_healpix_info(self) -> dict:
        """Get HEALPix-specific information.

        Returns
        -------
        info : dict
            Keys: ``nside``, ``npix_total``, ``pixel_area_sr``,
            ``pixel_area_arcmin2``, ``resolution_arcmin``,
            ``resolution_deg``, ``max_pixel_radius_deg``.

        """
        return {
            "nside": self.nside,
            "npix_total": self.hp.nside2npix(self.nside),
            "pixel_area_sr": self.hp.nside2pixarea(self.nside),
            "pixel_area_arcmin2": (
                self.hp.nside2pixarea(self.nside, degrees=True) * 3600
            ),
            "resolution_arcmin": self.hp.nside2resol(self.nside, arcmin=True),
            "resolution_deg": self.actual_angular_resolution,
            "max_pixel_radius_deg": np.rad2deg(self.hp.max_pixrad(self.nside)),
        }

class GeodesicBuilder(BaseGridBuilder):
    """Geodesic grid based on a subdivided icosahedron.

    The sphere is tessellated into triangular cells by starting with an
    icosahedron (20 equilateral triangles) and recursively subdividing each
    triangle into four smaller triangles.  All vertices are projected back
    onto the unit sphere after each subdivision step, so the final cells are
    *spherical* triangles.  The grid has no polar singularity and provides
    near-uniform cell areas, though strict equal-area is *not* guaranteed —
    cell areas vary by a few percent depending on how they inherit the
    icosahedral symmetry axes.

    Coordinate convention (physics / GNSS)
    ---------------------------------------
    * phi  ∈ [0, 2π)  – azimuthal angle from North, clockwise (navigation convention)
    * theta ∈ [0, π/2] – polar angle from zenith (0 = straight up,
      π/2 = horizon)

    Cell centres are computed as the 3D Cartesian mean of the three vertices,
    re-normalised onto the unit sphere.

    What ``angular_resolution`` means
    ----------------------------------
    ``angular_resolution`` is **not** used directly as a cell size.  Instead it
    is used only when ``subdivision_level`` is *not* explicitly supplied, to
    *estimate* an appropriate subdivision level.  The heuristic targets an
    approximate triangle edge length of ``2 × angular_resolution``::

        target_edge ≈ 2 × angular_resolution   (degrees)
        subdivision_level = ceil(log₂(63.4 / target_edge))

    The number 63.4° is the edge length of a regular icosahedron inscribed in
    a unit sphere.  Each subdivision halves the edge length, so the actual
    edge length at level *n* is approximately::

        edge ≈ 63.4° / 2ⁿ   (degrees)

    The total number of triangles on the **full sphere** is ``20 × 4ⁿ``.
    Roughly half fall in the northern hemisphere (exact count depends on
    the hemisphere boundary).

    Mathematical construction
    -------------------------
    1. **Icosahedron** – 12 vertices placed at the intersections of three
       mutually perpendicular golden-ratio rectangles, normalised to the
       unit sphere.  20 triangular faces connect them.
    2. **Subdivision** – each triangle is split into 4 by inserting edge
       midpoints.  Each midpoint is projected onto the unit sphere
       (re-normalised) before the next subdivision.  This is repeated
       ``subdivision_level`` times.
    3. **Hemisphere filter** – faces are kept if *any* of their three
       vertices satisfies ``theta ≤ π/2 − cutoff_theta``.  Consequently,
       boundary triangles that straddle the horizon *are* included and
       extend slightly below it.
    4. **Phi wrapping** – for triangles that straddle the 0/2π azimuthal
       boundary, vertex phis below π are shifted by +2π before computing
       bounding-box limits, then wrapped back.

    Parameters
    ----------
    angular_resolution : float
        Approximate angular resolution in degrees.  Used only to derive
        ``subdivision_level`` when that parameter is not given explicitly.
    cutoff_theta : float
        Elevation mask angle in degrees.  Triangles are excluded only if
        *all* their vertices are below this elevation.
    subdivision_level : int or None
        Number of icosahedral subdivisions.  If ``None``, estimated from
        ``angular_resolution``.  Typical range 0–5.
    phi_rotation : float
        Rigid azimuthal rotation applied after construction, in degrees.

    Notes
    -----
    The ``theta_lims``, ``phi_lims``, and ``cell_ids`` fields of the returned
    ``GridData`` are *synthetic* evenly-spaced arrays kept only for interface
    compatibility with ring-based grids.  They do **not** describe the actual
    triangular cell layout.  Use the ``geodesic_vertices`` column and the
    ``vertices`` array in ``GridData.vertices`` for the true geometry.

    """

    def __init__(
        self,
        angular_resolution: float = 2,
        cutoff_theta: float = 0,
        subdivision_level: int | None = None,
        phi_rotation: float = 0,
    ) -> None:
        """Initialize the geodesic grid builder.

        Parameters
        ----------
        angular_resolution : float, default 2
            Angular resolution in degrees.
        cutoff_theta : float, default 0
            Maximum polar angle cutoff in degrees.
        subdivision_level : int | None, optional
            Subdivision level override.
        phi_rotation : float, default 0
            Rotation angle in degrees.

        """
        super().__init__(angular_resolution, cutoff_theta, phi_rotation)
        self._triangles: np.ndarray | None = None

        if subdivision_level is None:
            target_edge_deg = angular_resolution * 2
            self.subdivision_level = max(
                0,
                int(np.ceil(np.log2(63.4 / target_edge_deg))),
            )
        else:
            self.subdivision_level = subdivision_level

        self._logger.info(
            f"Geodesic: subdivision_level={self.subdivision_level}, "
            f"~{20 * 4**self.subdivision_level} triangles"
        )

    def get_triangles(self) -> np.ndarray | None:
        """Return triangle vertex coordinates for visualization.

        Returns
        -------
        triangles : np.ndarray or None
            Array of shape ``(n_faces, 3, 3)`` where ``triangles[i]`` contains
            the three 3D unit-sphere vertices of triangle *i*.  ``None`` if
            the grid has not been built yet.

        """
        return self._triangles

    def get_grid_type(self) -> str:
        """Return the grid-type identifier string.

        Returns
        -------
        str
            ``"geodesic"``

        """
        return GridType.GEODESIC.value

    def _extract_triangle_vertices(
        self, vertices: np.ndarray, faces: np.ndarray
    ) -> np.ndarray:
        """Extract 3D vertex coordinates for each face.

        Parameters
        ----------
        vertices : np.ndarray
            All sphere vertices, shape ``(n_vertices, 3)``.
        faces : np.ndarray
            Face index array, shape ``(n_faces, 3)``.

        Returns
        -------
        triangles : np.ndarray
            Shape ``(n_faces, 3, 3)`` – three 3D vertices per face.

        """
        # Vectorized: use NumPy advanced indexing instead of loop
        return vertices[faces]

    def _build_grid(
        self,
    ) -> tuple[
        pl.DataFrame, np.ndarray, list[np.ndarray], list[np.ndarray], dict[str, Any]
    ]:
        """Build geodesic grid from subdivided icosahedron.

        Returns
        -------
        grid : pl.DataFrame
            One row per triangular cell.  Columns include phi, theta (centre),
            bounding-box limits, ``geodesic_vertices`` (3 vertex indices into
            the ``vertices`` array), and ``geodesic_subdivision``.
        theta_lims : np.ndarray
            Synthetic evenly-spaced theta limits (interface compatibility only).
        phi_lims : list[np.ndarray]
            Synthetic evenly-spaced phi limits (interface compatibility only).
        cell_ids : list[np.ndarray]
            Single-element list containing all cell ids.
        extra_kwargs : dict
            Contains ``vertices`` (shape ``(n_vertices, 3)``),
            ``vertex_phi``, and ``vertex_theta`` arrays for the full
            subdivided icosahedron.

        """
        vertices, faces = self._create_icosahedron()

        # Subdivide
        for _ in range(self.subdivision_level):
            vertices, faces = self._subdivide_mesh(vertices, faces)

        # Project to unit sphere
        vertices = vertices / np.linalg.norm(vertices, axis=1, keepdims=True)

        # Convert to spherical
        x, y, z = vertices[:, 0], vertices[:, 1], vertices[:, 2]
        theta = np.arccos(np.clip(z, -1, 1))
        phi = np.arctan2(y, x)
        phi = np.mod(phi, 2 * np.pi)

        # Filter to northern hemisphere
        hemisphere_mask = theta <= (np.pi / 2 - self.cutoff_theta_rad)

        # Filter faces
        valid_faces = []
        for face in faces:
            if any(hemisphere_mask[v] for v in face):
                valid_faces.append(face)

        if len(valid_faces) == 0:
            raise ValueError("No valid faces in hemisphere")

        valid_faces = np.array(valid_faces)

        # Create cells
        cells = []
        for face in valid_faces:
            v_indices = face
            face_phi = phi[v_indices]
            face_theta = theta[v_indices]

            # Handle phi wrapping for triangles crossing 0/2π boundary
            phi_range = np.ptp(face_phi)
            if phi_range > np.pi:
                # Triangle crosses the wraparound - unwrap relative to median
                ref_phi = np.median(face_phi)
                face_phi_unwrapped = face_phi.copy()
                # Unwrap angles that are > π away from reference
                mask_low = (ref_phi - face_phi_unwrapped) > np.pi
                mask_high = (face_phi_unwrapped - ref_phi) > np.pi
                face_phi_unwrapped[mask_low] += 2 * np.pi
                face_phi_unwrapped[mask_high] -= 2 * np.pi
                phi_min = float(np.min(face_phi_unwrapped) % (2 * np.pi))
                phi_max = float(np.max(face_phi_unwrapped) % (2 * np.pi))
            else:
                phi_min = float(np.min(face_phi))
                phi_max = float(np.max(face_phi))

            # Cell center - 3D Cartesian mean
            face_vertices_3d = vertices[v_indices]
            center_3d = np.mean(face_vertices_3d, axis=0)
            center_3d = center_3d / np.linalg.norm(center_3d)

            center_theta = np.arccos(np.clip(center_3d[2], -1, 1))
            center_phi = np.arctan2(center_3d[1], center_3d[0])
            center_phi = np.mod(center_phi, 2 * np.pi)

            # Cell bounds (theta from vertices, phi already computed above)
            theta_min = float(np.min(face_theta))
            theta_max = float(np.max(face_theta))

            cells.append(
                {
                    "phi": center_phi,
                    "theta": center_theta,
                    "phi_min": phi_min,
                    "phi_max": phi_max,
                    "theta_min": theta_min,
                    "theta_max": theta_max,
                    "geodesic_vertices": v_indices.tolist(),
                    "geodesic_subdivision": self.subdivision_level,
                }
            )

        grid = pl.DataFrame(cells).with_columns(
            pl.int_range(0, pl.len()).alias("cell_id")
        )

        extra_kwargs: dict[str, Any] = {
            "vertices": vertices,
            "vertex_phi": phi,
            "vertex_theta": theta,
        }

        theta_lims = np.linspace(0, np.pi / 2, 10)
        phi_lims = [np.linspace(0, 2 * np.pi, 20) for _ in range(len(theta_lims))]
        cell_ids_list = [np.arange(grid.height)]

        self._triangles = self._extract_triangle_vertices(vertices, faces)

        return grid, theta_lims, phi_lims, cell_ids_list, extra_kwargs

    def _create_icosahedron(self) -> tuple[np.ndarray, np.ndarray]:
        """Create a unit-sphere icosahedron.

        Returns
        -------
        vertices : np.ndarray
            Shape ``(12, 3)`` – vertices on the unit sphere.
        faces : np.ndarray
            Shape ``(20, 3)`` – integer vertex indices per triangular face.

        """
        phi_golden = (1 + np.sqrt(5)) / 2

        vertices = np.array(
            [
                [-1, phi_golden, 0],
                [1, phi_golden, 0],
                [-1, -phi_golden, 0],
                [1, -phi_golden, 0],
                [0, -1, phi_golden],
                [0, 1, phi_golden],
                [0, -1, -phi_golden],
                [0, 1, -phi_golden],
                [phi_golden, 0, -1],
                [phi_golden, 0, 1],
                [-phi_golden, 0, -1],
                [-phi_golden, 0, 1],
            ],
            dtype=np.float64,
        )

        vertices = vertices / np.linalg.norm(vertices, axis=1, keepdims=True)

        faces = np.array(
            [
                [0, 11, 5],
                [0, 5, 1],
                [0, 1, 7],
                [0, 7, 10],
                [0, 10, 11],
                [1, 5, 9],
                [5, 11, 4],
                [11, 10, 2],
                [10, 7, 6],
                [7, 1, 8],
                [3, 9, 4],
                [3, 4, 2],
                [3, 2, 6],
                [3, 6, 8],
                [3, 8, 9],
                [4, 9, 5],
                [2, 4, 11],
                [6, 2, 10],
                [8, 6, 7],
                [9, 8, 1],
            ],
            dtype=np.int64,
        )

        return vertices, faces

    def _subdivide_mesh(
        self, vertices: np.ndarray, faces: np.ndarray
    ) -> tuple[np.ndarray, np.ndarray]:
        """Subdivide each triangle into 4 smaller triangles.

        Each edge midpoint is computed, normalised onto the unit sphere, and
        cached so that shared edges produce only one new vertex.

        Parameters
        ----------
        vertices : np.ndarray
            Current vertex array, shape ``(n_vertices, 3)``.
        faces : np.ndarray
            Current face array, shape ``(n_faces, 3)``.

        Returns
        -------
        new_vertices : np.ndarray
            Expanded vertex array, shape ``(n_vertices + n_new_midpoints, 3)``.
        new_faces : np.ndarray
            New face array, shape ``(4 × n_faces, 3)``.

        """
        new_faces = []
        edge_midpoints: dict[tuple[int, int], int] = {}

        def get_midpoint(v1: int, v2: int) -> int:
            """Return midpoint vertex index for an edge.

            Parameters
            ----------
            v1 : int
                First vertex index.
            v2 : int
                Second vertex index.

            Returns
            -------
            int
                Index of the midpoint vertex.

            """
            edge = tuple(sorted([v1, v2]))
            if edge not in edge_midpoints:
                edge_midpoints[edge] = len(vertices) + len(edge_midpoints)
            return edge_midpoints[edge]

        for face in faces:
            v0, v1, v2 = face

            m01 = get_midpoint(v0, v1)
            m12 = get_midpoint(v1, v2)
            m20 = get_midpoint(v2, v0)

            new_faces.extend(
                [
                    [v0, m01, m20],
                    [v1, m12, m01],
                    [v2, m20, m12],
                    [m01, m12, m20],
                ]
            )

        n_original = len(vertices)
        n_new = len(edge_midpoints)
        final_vertices = np.zeros((n_original + n_new, 3))
        final_vertices[:n_original] = vertices

        for edge, idx in edge_midpoints.items():
            v1, v2 = edge
            midpoint = (vertices[v1] + vertices[v2]) / 2
            midpoint = midpoint / np.linalg.norm(midpoint)
            final_vertices[idx] = midpoint

        return final_vertices, np.array(new_faces)